Equilibrium swelling and universal ratios in dilute polymer solutions: Exact Brownian 
dynamics simulations for a delta function excluded volume potential 
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A narrow Gaussian excluded volume potential, which tends to a J-function repulsive potential in 
the limit of a width parameter d* going to zero, has been used to examine the universal consequences 
of excluded volume interactions on the equilibrium and linear viscoelastic properties of dilute poly- 
mer solutions. Brownian dynamics simulations data, acquired for chains of finite length, has been 
extrapolated to the limit of infinite chain length to obtain model independent predictions. The 
success of the method in predicting well known aspects of static solution properties suggests that 
it can be used as a systematic means by which the influence of solvent quality on both equilibrium 
and non-equilibrium properties can be studied. 

PACS numbers: 83.10.Mj, 61.25.Hq, 83.60.Bc, 83.80.Rs 
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There is a growing recognition of the importance of sol- 
vent quality in determining the rheological properties of 
dilute polymer solutions. Though there are relatively few 
systematic experimental investigations of the influence of 
solvent quality, nevertheless there is sufficient evidence 
that material functions, in both shear and extensional 
flows, are significantly different from each other in good 
and 6>-solvents 0, 0> Q . While the theoretical descrip- 
tion of the influence of solvent quality on the equilibrium 
behavior of dilute polymer solutions has been a major 
area of research [J, the description of its influence on 
non- equilibrium behavior is still in its infancy. Most cur- 
rent non-equilibrium theories do not attempt to develop 
a unified framework that is also applicable at equilib- 
rium, and the wealth of experience and insight that has 
been gained so far in the development of static theories 
is rarely used in the description of dynamic behavior. 
In this paper, we introduce an approach that provides a 
consistent means of describing both regimes of behavior, 
and as a first step towards establishing the usefulness of 
the methodology, we demonstrate that the method is ca- 
pable of reproducing well known universal static results. 
We also show that, even in the case of well established 
results, the method is capable of providing new and ad- 
ditional insights. Furthermore, as an example of the ver- 
satility of the approach, the influence of solvent quality 
on a universal ratio of linear viscoelastic properties, is 
reported here for the first time. 

One of the most important results of the experimental 
investigation of the static behavior of dilute polymer so- 
lutions, has been the discovery that various properties — 
both in ^-solvents and in good solvents — exhibit power 
law behavior when the molecular weight of the dissolved 
polymer is sufficiently large. It is perhaps less com- 
monly known that even when the molecular weight is not 
very large, a type of scaling still persists, which enables 
the description of behavior with the help of a single pa- 
rameter. An illustrative example is the behavior of the 
root mean square radius of gyration, R g . Light scat- 



tering measurements indicate that for molecules of suffi- 
ciently large molecular weight M, R g obeys the power law 
i?g = aM v , where the pre-factor a depends on the par- 
ticular polymer-solvent system, but the exponent v does 
not. In a ^-solvent, the universal exponent has a value, 
v = 0.5, while in a good solvent, v = 0.592 ±0.003 0. A 
remarkable experimental observation is that, away from 
the large molecular weight asymptotic limit, the depen- 
dence of the mean size of the polymer molecule on both 
the temperature, T, of the solution and M , can be com- 
bined into a single variable, r(T)\/M, where, r(T) is a 
simple function of temperature, r = (1 — Tg/T), with 
Tq denoting the 0-temperature. A vast amount of equi- 
librium data for a variety of polymer solvent systems re- 
veals that, when the swelling cv g of a polymer molecule — 
defined as the ratio of the radius of gyration in a good 
solvent to the radius of gyration in a 0-solvent — is plotted 
as a function of ry/M (with each polymer-solvent system 
being shifted horizontally by a suitable constant factor), 
universal behavior is displayed, for all values of the molec- 
ular weight M, and temperature T > T e [EIS0- 

The entire range of behavior exhibited by static solu- 
tion properties has been successfully predicted by apply- 
ing renormalization group (RG) methods [4j. In these 
theories, the macromolecule is frequently represented by 
a coarse-grained model, such as a bead-spring chain, con- 
sisting of N beads connected together by (N— 1) Hookean 
springs (with a spring constant H). The presence of ex- 
cluded volume interactions is then taken into account by 
assuming the existence of a Dirac delta function repul- 
sive potential, that acts pairwise between the beads of 
the chain 8], 



E{r^)=v(T) k B TS{r^) 



(1) 



where, v(T) is the excluded volume parameter, fee is 
Boltzmann's constant, and r utl is the vector connecting 
beads v and fi. In the limit T — ► Tg, it can be shown that 
v depends on T through the relation, v(T) = vqt(T), 
where vq is a constant ||. If a length scale, £ = 
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y/knT/H, and a non-dimensional strength of excluded 

volume interactions, z* = v (2tt£ 2 ) 3 ^ 2 are defined, then 
both the temperature and the chain length dependence 
can be combined into the single non-dimensional parame- 
ter, z = z*yN. Completely consistent with experimental 
observations (since z oc t\/M), RG theories predict the 
existence of power laws in the limit z — > (correspond- 
ing to (9-solvents), and the excluded volume limit z — > oo 
(corresponding to good solvents), and the existence of 
scaling functions, that depend only on the parameter z, 
which accurately describe the crossover behavior between 
these two asymptotic limits 4]. 

In spite of the success of RG theories, it must be borne 
in mind that they are approximate theories, dependent 
to some extent on the order of the original perturbation 
calculation on which they are based. Exact numerical re- 
sults for the excluded volume problem can, on the other 
hand, be obtained by either Monte Carlo simulations of 
self avoiding walks on lattices 0, Eli or by off-lattice 
Monte Carlo simulations Since S- functions can- 

not be used in numerical investigations, off-lattice Monte 
Carlo simulations are typically based on excluded volume 
potentials with a finite range of excluded volume inter- 
actions, such as the Lennard- Jones potential. An impor- 
tant point to be noted, however, is that the parameter 
z, which is the true measure of solvent quality, does not 
appear naturally when potentials like the Lennard-Jones 
potentials are used. Instead, simulation results are func- 
tions of the parameters that characterize the potential. 
It is possible, however, to infer the value of z from sim- 
ulation results by carrying out a data shifting procedure 
similar to that used to bring experimental data onto a 
universal curve, i.e, data for various values of chain length 
JV, and various distances from the ^-temperature, can be 
brought to lie on a master curve 11]. The procedure is 
quite involved and even the ^-temperature must first be 
estimated from simulations. Furthermore, Monte Carlo 
methods cannot be used to estimate rheological prop- 
erties in general flow fields, and so it is not clear how 
this procedure can be extended to systematically exam- 
ine the influence of solvent quality on non-equilibrium 
properties. 

The self-similar character of polymer molecules, which 
is responsible for the universal behavior exhibited by 
polymer solutions at equilibrium, is also responsible 
for the universal behavior displayed away from equilib- 
rium ||. An important and fundamental challenge for 
molecular theories is to verify the existence of this uni- 
versal behavior. The universal consequences of excluded 
volume effects away from equilibrium has been exam- 
ined by Ottinger and co-workers fl^ . fl3| with the help of 
RG methods. Since these methods lead to approximate 
predictions, it is important to be able to assess their ac- 
curacy. Further, these papers were focussed on predict- 
ing the universal behavior in the excluded volume limit, 



and cross-over scaling functions for rheological properties 
were not reported. 

In this paper, we introduce a procedure by which uni- 
versal cross-over scaling functions and asymptotic behav- 
ior in the excluded volume limit can be obtained for a 
(5-function excluded volume potential, both at equilib- 
rium and away from it. In other words, a means has 
been obtained by which an exact numerical solution to 
the Edwards model, which is of fundamental importance 
to polymer physics, can be found. Using the suggested 
methodology, provided sufficient computational time is 
used, properties predicted by the Edwards model can be 
found to a high degree of accuracy. As a result, the in- 
fluence of solvent quality, measured directly in terms of 
the parameter z, can be systematically examined, and 
results obtained by approximate RG calculations may be 
verified. 

A detailed discussion of the theoretical framework used 
here may be found in earlier papers by Prakash |l4i 1 1 ol| . 
Essentially, macroscopic properties are obtained by ex- 
ploiting the mathematical equivalence of diffusion equa- 
tions in polymer configuration space and stochastic dif- 
ferential equations for the polymer configuration. Since 
averages calculated from stochastic trajectories are iden- 
tical to averages calculated from distribution functions, 
stochastic trajectories are generated by numerically in- 
tegrating the appropriate stochastic differential equa- 
tion with the help of Brownian dynamics simulations 
(BDS) 0. 

A key ingredient in our approach is the use of a narrow 
Gaussian potential (NGP), 

^> = (^)%t„p{-^^} ( 2) 

where, al* is an additional non-dimensional parameter 
that measures the range of excluded volume interaction. 
The NGP is a means of regularizing the (5-function po- 
tential since it reduces to the (^-potential in the limit of 
d* — > 0. The central hypothesis behind the use of this 
potential is that when data accumulated for chains of fi- 
nite length is extrapolated to the infinite chain length 
limit, universal consequences of the presence of excluded 
volume interactions will be revealed, independent of the 
value of d* . Initial attempts to carry out this program 
were unsuccessful because of the computational intensity 
of BDS. The essential soundness of the procedure has, 
however, been established by obtaining universal predic- 
tions, both at equilibrium and in the linear viscoelastic 
limit 1171 . and at finite shear rates in steady simple shear 
flow [l£j , with an approximate version of the theory — the 
so called Gaussian approximation. The parameter d* has 
been shown to be irrelevant in the long chain limit, since 
it always appears in the theory as the ratio d* /\f~N . 

An exact numerical solution of the model has been 
obtained in this work, and the computational difficulties 
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encountered earlier have been overcome, by adopting two 
new strategies (apart from the obvious one of paralleliz- 
ing the code). First, a novel trajectory in the (d*,z*) 
parameter space has been used to reach the asymptotic 
limit. Second, a variance reduction procedure based on 
the method of control variates ^(| , with the Rouse model 
as the control variate, has been implemented. While 
these strategies, which are discussed in greater detail else- 
where do not alter the scaling of CPU time with N 
(which continues to increase roughly as A 3 ), they ensure 
that the noise is sufficiently small to permit an accurate 
extrapolation of finite chain data, even for fairly small 
chains. 

Two different approaches have been adopted here to 
obtain results in the limit N — > oo. In the first, z is kept 
constant while simulations are carried out with increas- 
ing values of N. Extrapolation of the accumulated data 
to N — > oo implies the simultaneous limit, z* — > (or 
equivalently, T — ► Tg). As a result, crossover functions 
for systems close to the ^-temperature are obtained. In 
the second approach, z* is maintained constant (which 
is equivalent to T being constant) while N is increased 
systematically. Extrapolation to N — > oo in this case 
leads to the prediction of behavior in the excluded vol- 
ume limit, z — > oo. 

The crossover scaling function for the swelling a^, as 
a function of z, is displayed in Fig. ^a). The exper- 
imental data of Miyaki and Fujita [6j for polystyrene 
of various molecular weights and temperatures has been 
shifted to lie on the simulation data, using a single con- 
stant horizontal shift factor for each of the two sol- 
vents 18] . The remarkable agreement of the predicted 
equilibrium swelling with experimental data, similar to 
that obtained earlier by RG theory 0, has been ob- 
tained by using bead-spring chains with a maximum of 
just 36 beads. The BDS data has been curve-fitted with 
an equation commonly used to fit results of RG analy- 
sis 0, c*g = (l + az + bz 2 +cz 3 ) 14 "" 2 " 3 , with the coef- 
ficient a kept fixed at its first order perturbation value . 
It leads to the prediction v = 0.6004 ± 0.0004, which is 
about 2% greater than the value v — 0.588 predicted by 
RG theory. The maximum difference between the com- 
puted data and the curve fit was less than 0.6%. As 
shown below, we find that v estimated in the excluded 
volume limit, with the present procedure, is also about 
2% more than 0.588. 

The efficacy of the present approach is demonstrated 
in Fig.^b), where the crossover behavior of the univer- 
sal ratio constructed from the zero shear rate first nor- 
mal stress difference coefficient and the zero shear rate 
viscosity, Uqj v — n p kBT^> i,o/Vp,o (with n p denoting the 
number density of polymers), is displayed. To our knowl- 
edge, this crossover behavior has not been reported else- 
where before. The only existing universal prediction is a 
RG calculation in the excluded volume limit (indicated 
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FIG. 1: Crossover scaling functions, (a) Swelling of the radius 
of gyration, a|, versus the strength of excluded volume inter- 
action, z = z*(T)\ /r N. Experimental data for polystyrene in 
the two solvents is that of Miyaki and Fujita 0- (b) Universal 
ratio U<sr) as a function of z. The solid line is the prediction 
in the excluded volume limit, while the dashed line is the 
prediction of RG theory. 

in the figure by the dashed line), obtained by refining a 
first order perturbation expansion 

In the Rouse model, = constant = 0.8. It is also 
expected to be constant, independent of the polymer- 
solvent system, in the excluded volume limit. The figure 
suggests that the crossover function approaches the value 
predicted by RG theory, U^ n = 0.6288, as z — > oo. On 
the other hand, as indicated by the solid line, a signifi- 
cantly different value of U*^ = 0.77 ± 0.01 is predicted 
by the present approach in the excluded volume limit. It 
is of interest to see if a RG calculation based on a higher 
order perturbation theory would lead to a prediction in 
agreement with the present prediction. A similar differ- 
ence, between the large z limit suggested by the crossover 
behavior, and the value obtained in the excluded volume 
limit, has been observed previously in the case of the uni- 
versal ratio Ur — a?/ a?, where a 2 is the swelling of the 
end-to-end vector. The behavior of Ur predicted by the 
present approach will be discussed elsewhere 0] . 

Although the scaling of R g with chain length N obeys 
a power law only in the limit of long chains, one can 
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FIG. 2: The exponent v e s, in the expression R s ~ N" ef! , as a 
function of 2, and in the inset, as a function 2*. 



make the ansatz, i? g ~ N" et! , with ;/ e ff representing an 
effective exponent that approaches its critical value as 
N — ► 00. By carrying out high precision Monte Carlo 
studies of lattice self avoiding walks, Sokal and cowork- 
ers have observed that is c g approaches its asymptotic 
value from above, in contrast to modern two-parameter 
theories for the continuum Edwards model that predict 
that the asymptotic value is reached from below. This 
has prompted them to suggest that these theories do 
not correspond to good solvents in general, but instead 
to an infinitesimal region just above the theta tempera- 
ture. Schafer and co-workers 0, [ljj have shown, by de- 
veloping an alternative renormalization group method, 
and by Monte Carlo simulations, that both approaches 
to the asymptotic limit are not incompatible due to 
the existence of a two branched structure to the so- 
lution of the excluded volume problem. Their results 
indicate that the excluded volume limit is approached 
on two distinct branches, depending on the whether 
z* < Zf (the weak-coupling branch), or z* > z\ (the 
strong-coupling branch), where z\ denotes the fixed point 
for the parameter z* . While the existence of a dual 
branched structure has only recently been elucidated, 
the various curves in Fig. [21 for v e g as a function of 
z, show that Brownian dynamics simulations readily re- 
veal the presence of this structure. Indeed, as will be 
discussed below, a unique way can be developed by 
which the distinctive structure of the solution may be 
exploited to find both the fixed point and the critical 
exponent. The curves in Fig. [3 have been obtained (i) 



by fitting or versus z data, for z* < 



(l + az + 5z 2 ) 2 " ; (ii) by fitting the data for z* > z\ 
with a\ = az iu - 2 [l -&2-H (l + c/z-l 1+m l)] ; and (iii) 
by using the fact that v e s = 0.5 + 0.25(91nag/<9mz). 
In both the curve-fitting expressions used for parameter- 
izing the data (whose forms have been suggested earlier [17 
by Grassberger et al. 0), the value of v has been set E l8 



4 f , with cVg 



equal to 0.6. 

On the weak-coupling branch (strong-coupling 
branch), v e g first increases (decreases) rapidly before 
approaching the asymptotic value very gradually. We 
consequently expect a plot of v c r vs z* , for constant 
values of 0?* and z, to show a point of inflection at the 
fixed point. This is clearly seen in the inset of Fig. [21 
which describes the change in v c g as a function of z*, 
for a range of values of z* spanning the fixed point. 
The figure suggests that the fixed point lies in the range 
0.28 < z\ < 0.3, and that the value of v e s in this range 
is 0.6004 ± 0.0002, which is expected to be close to the 
critical value, since asymptotic behavior is attained at 
relativelysmall values of TV as z* — > z\. Schafer and co- 
workers li-lni have estimated z\ = 0.364, while recently 
Prakash [l5j has shown in the Gaussian approximation 
that, Zf « 0.35. It is an exciting prospect to see if this 
two branched structure persists into the non-equilibrium 
regime — an issue that is easily examinable within the 
framework introduced in this work. 
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